AFML-TR-67-121 
PART  VII 

frVOlU'BlZ 


EVALUATION  OF  MOLECULAR  WEIGHT 
FROM  EQUILIBRIUM  SEDIMENTATION 


PART  VII  MATHEMATICAL  ANALYSIS  OF  THE  REGULARIZATION 
TECHNIQUE  INCORPORATED  INTO  QUADRATIC  PROGRAMING 


DONALD  R.  WIFF 
MATATIAHU  T.  GEHATIA 
THOMAS  E.  DUVALL 


TECHNICAL  REPORT  AFML-TR-67-121,  PART  VII 


DECEMBER  1972 


,  7/1* 

Approved  for  public  release;  distribution  unlimited. 

f  Jr  / 


V  / 


{y/  /.£> ; 


AIR  FORCE  MATERIALS  LABORATORY 
AIR  FORCE  SYSTEMS  COMMAND 
WRIGHT-PATTERSON  AIR  FORCE  BASE,  OHIO 

BEST  AVAILABLE  COPY 

3jDo^o3aa  m 


NOTICE 


When  Government  drawings,  specifications,  or  other  data  are  used  for  any  purpose 
other  than  in  connection  with  a  definitely  related  Government  procurement  operation, 
the  United  States  Government  thereby  incurs  no  responsibility  nor  any  obligation 
whatsoever;  and  the  fact  that  the  government  may  have  formulated,  furnished,  or  in 
any  way  supplied  the  said  drawings,  specifications,  or  other  data,  is  not  to  be  regarded 
by  implication  or  otherwise  as  in  any  manner  licensing  the  holder  or  any  other  person 
or  corporation,  or  conveying  any  rights  or  permission  to  manufacture,  use,  or  sell  any 
patented  invention  that  may  in  any  way  be  related  thereto. 


Copies  of  this  report  should  not  be  returned  unless  return  is  required  by  security 
considerations,  contractual  obligations,  or  notice  on  a  specific  document, 

AIR  FORCE/567 80/22  February  1973—  100 


AFML-TR-67-121 
PART  VII 


EVALUATION  OF  MOLECULAR  WEIGHT 
FROM  EQUILIBRIUM  SEDIMENTATION 

PART  VII  MATHEMATICAL  ANALYSIS  OF  THE  REGULARIZATION 
TECHNIQUE  INCORPORATED  INTO  QUADRATIC  PROGRAMING 


DONALD  R.  WIFF 
MATATIAHU  T.  GEHATIA 
THOMAS  E.  DUVALL 


Approved  for  public  release;  distribution  unlimited. 


FOREWORD 


AFML-TR-67-1 21 
PART  VII 


This  report  was  prepared  by  the  Polymer  Branch  of  the  Nonmetallic 
Materials  Division.  The  work  was  initiated  under  Project  No.  7342, 
"Fundamental  Research  on  Macromolecular  Materials  and  Lubrication 
Phenomena/'  Task  No.  734203,  "Fundamental  Principles  Determining  the 
Behavior  of  Macromolecules/1  Subtask  No.  734203-0,5,  "Physical  Chemistry 
of  High  Polymers",  with  Dr.  M.  T.  Gehatia  acting  as  subtask  scientist. 
Coauthors  are  Mr.  T.  E.  Duvall,  ASD  Computer  Science  Center  (4950/VNCS), 
and  Dr.  D.  R.  Wiff,  Research  Institute,  University  of  Dayton,  The  work 
was  administered  under  the  direction  of  the  Air  Force  Materials 
Laboratory,  Air  Force  Systems  Command,  Wright-Patterson  Air  Force  Base, 
Ohio. 

The  report  covers  research  conducted  from  September  1971  to  May  1972. 
The  manuscript  was  released  by  the  authors  in  June  1972  for  publication 
as  a  technical  report. 

This  technical  report  has  been  reviewed  and  is  approved. 


R.  L.  VAN  DEUSEN 
Chief,  Polymer  Branch 
Nonmetallic  Materials  Division 
Air  Force  Materials  Laboratory 


ii 


ABSTRACT 


The  equation  relating  molecular  weight  distribution  of  a  polymer  to 
the  experimental  function  of  concentration  appearing  in  equilibrium 
sedimentation  with  the  ultracentrifuge  is  nonsolvable  because  it  is  an 
Improperly  Posed  Problem  in  the  Hadamard  sense.  For  a  simple  distri¬ 
bution  this  equation  has  been  solved  by  applying  a  method  of  regular¬ 
ization.  To  solve  a  nonsymmetrical  bimodal  and  a  trimodal  distribution, 
the  technique  of  regularization  had  to  be  incorporated  into  a  linear 
programming.  In  the  current  work  the  regularization  technique  has  been 
incorporated  into  quadratic  programming.  This  new  combined  method  proved 
to  be  more  adequate  to  solve,  also  more  complex  distributions  such  as 
tri-,  tetra-,  and  pentamodal .  In  addition  this  technique  is  cheaper, 
because  it  requires  less  computer  time  than  the  regularization  incor¬ 
porated  into  linear  programming. 
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SECTION  I 
INTRODUCTION 

The  increased  use  of  polymeric  materials  by  the  U.S.  Air  Force  has 
placed  an  ever  increasing  demand  upon  the  reliability  e.  g.,  strength,  of 
these  materials.  Many  bulk  property  characteristics:  density,  shear 
modulus,  stress  modulus,  high  temperature  resistance,  tenacity,  etc.  are 
dependent  upon  the  distribution  of  molecular  weights  of  the  macromolecular 
chains  composing  the  material. 

There  are  many  interacting  morphological  patterns  -  tie  molecules, 
degrees  of  crystallinity,  varying  degrees  of  order-which  manifest  various 
bonding  energies  and/or  intra-molecular  interactions.  These  affect  the 
strength  of  a  polymeric  material  in  the  bulk  state.  However,  if  the 
molecular  weight  is  too  low,  the  strength  can  be  affected  as  a  result  of 
pure  thermal  (Brownian)  motion.  An  extremely  high  molecular  weight 
might,  on  the  other  hand,  inhibit  relaxation  or  even  hinder  the  process- 
ability  of  the  material.  If  molecular  weight  affects  the  final  bulk 
state  properties  to  such  a  degree,  a  distribution  of  molecular  weights 
adds  another  variable  that  can  greatly  affect  the  reliability  of  these 
materials. 

For  these  reasons,  a  mathematical  procedure  for  obtaining  a  molecular 
weight  distribution  (MWD)  from  equilibrium  sedimentation  data  was 
necessary. 

There  exist  differential  and  integral  equations  describing  important 
physical  or  technological  systems  which  in  general  cannot  be  solved  by 
usual  mathematical  means  and  even  by  approximation  because  they  belong  to 
the  class  of  improperly  posed  problems  (IPP).  To  this  class  also  belongs 
the  equation  which  relates  MWD  to  the  concentration  function  provided  by 
the  technique  of  equilibrium  sedimentation. 

The  notion  of  an  IPP  (improperly  posed  problem  or  incorrectly 
formulated  problem)  goes  back  to  Hadamard  (Reference  1)  in  conjunction 
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with  the  Cauchy  problems  of  potential  and  a  number  of  inverse  problems 
for  differential  and  integral  equations.  In  the  recent  decade  IPP's  have 
been  intensively  investigated.  The  following  considerations  with  respect 
to  ill-posedness  of  a  mathematical  problem  and  the  ways  leading  to  their 
solution  is  based  on  the  ideas  of  Phillips  (Reference  2),  John  (References 
3,  4),  Lavrentiev  (Reference  5),  Tikhonov  (References  6-21),  Ivanov 
(References  22-26),  and  others  (References  27-42).  Among  the  class  of 
IPP  there  exists  a  subclass  of  regularizable  IPP  which  can  be  solved  by 
applying  a  method  of  regularization. 

To  the  subclass  of  regularizable  IPP  also  belong  the  equation 
mentioned  before  associated  with  MWD  determination  via  equilibrium  sedi¬ 
mentation.  Because  of  the  need  to  correlate  a  MWD  with  physical  and 
mechanical  properties  of  synthetic  polymers  an  attempt  has  been  made  in 
this  laboratory  to  solve  this  particular  equation.  The  progress  of  this 
work  has  been  described  in  a  series  of  technical  reports  AFML-TR-67-1 21 , 
Parts  I  through  VI.  The  first  attempts  to  derive  an  MWD  from  these 
equations  without  using  the  regularization  technique  were  unsatisfactory 
(Parts  I  through  HI).  In  part  IV  regularization  was  successfully  applied 
and  good  results  were  obtained  in  case  of  a  unimodal  distribution.  To 
solve  more  complex  distributions,  such  as  symmetrical  and  asymmetrical 
bimodal  and  symmetrical  trimodal ,  regularization  was  incorporated  into 
a  linear  programing  algorithm  (Part  V).  In  Part  VI  this  method  was 
experimentally  verified.  An  artificial  and  a  priori  known  distribution 
of  polystyrene  samples  was  investigated.  The  resulting  distribution  was 
in  very  good  agreement  with  the  one  artificially  prepared. 

This  regularization  -  linear  programing  technique  seemed  limited  to  a 
maximum  trimodal  multiplicity.  In  addition,  a  large  amount  of  valuable 
digital  computer  time  was  consumed  in  search  for  appropriate  regularizing 
parameters. 

Therefore  the  present  report  (Part  VII)  extends  the  previously 
discussed  modifications  to  include  regularization  into  quadratic 
programing.  The  required  computation  time  is  greatly  diminished  and  the 


2 


AFML-TR-67-1 21 
PART  VII 


multiplicity  capable  of  being  resolved  now  includes  a  tetramodal  MWD. 

This  paper  is  divided  into  sections,  such  that,  once  the  mathematical 
definition  of  an  ill -posed  problem  is,  specified,  the  technique  of 
regularization  used  later  in  the  discussion,  will  be  explained.  This 
technique  leads  to  better  results  if  it  is  incorporated  in  the  quadratic 
programing.  However,  before  discussing  this  latest  refined  combination 
of  methods,  a  brief  discussion  of  a  quadratic  programing  technique 
follows.  Then  the  actual  Fredholm  Integral  of  the  First  Kind  used  along 
with  examples  of  its  ill-posedness  is  illustrated.  Finally,  the  incor¬ 
poration  of  regularization  into  quadratic  programing  with  its  application 
to  a  specific  kernel  will  be  presented. 

The  preceding  AFML-TR-67-1 21  reports  previously  referred  are: 


Part 

I  , 

M. 

T. 

Part 

II  , 

M. 

T. 

Part 

HI, 

R. 

R. 

Part 

IV  , 

M. 

T. 

Part 

V  , 

D. 

R. 

Part 

VI  , 

M. 

T. 

Gehatia  (June  1967). 

Gehatia  and  D.  R.  Wiff  (April  1969). 

Jurick,  D.  R.  Wiff  and  M.  T.  Gehatia  (May  1970). 
Gehatia  and  D.  R.  Wiff  (August  1970). 

Wiff  and  M.  T.  Gehatia  (February  1971). 

Gehatia  and  D.  R.  Wiff  (November  1971). 
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SECTION  II 

ILL-POSED  PROBLEM  AND  REGULARIZATION 

Let  F  and  U  be  some  complete  metric  spaces.  Let  Af  be  a  function 
with  domain  of  definition  F  and  the  range  of  values  U.  Consider  the 
equation 


Af  =  U  =  Q  [  £,  f(m)] 


The  problem  of  solving  Equation  1  for  a  set  {f}  given  a  set  {u }  and 
knowing  the  functional  form  of  A  is  a  properly  posed  problem  if  the 
following  conditions  are  satisfied: 

(la)  The  solution  of  Equation  1  exists  for  any  ueU. 

(lb)  The  solution  of  Equation  1  is  unique  in  F. 

(l c )  The  solution  of  Equation  1  depends  continuously  on  u  in  the 

metrics  of  F  and  U.  In  such  a  case  there  exists  a  function  Ou  defined 
and  continuous  over  all  U,  and  0  is  an  inverse  operator  of  A,  where 

Ou  =  A-1  u  =  f  =  R  [m,  u(f)  ]  (2) 


If  even  one  of  the  conditions  (la),  (lb)  or  (1c)  is  not  satisfied 
[u  =  Af]  is  an  IPP.  In  such  a  case  the  function  0  either  does  not  exist 
or  it  is  not  stable  and  not  reliable.  Many  expressions  of  mathematical 
physics  include  linear  operations.  In  this  case  U  and  F  are  Banach 
Spaces  and  A  is  a  linear  operator.  The  Banach  Spaces  U  and  F  encountered 
in  most  cases  are  the  known  functional  spaces  CA,  L  ,  Wp,  H^,  S  ,  .  .with 
the  carriers  in  some  n-dimensional  space  of  the  independent  variables  or 
on  any  part  of  th'e  spaces  of  independent  variables.  The  first  require¬ 
ment  of  correctness  is  that  the  problem  under  consideration  should  not 
be  overdetermined*,  second  that  the  solution  is  unique,  since  the  right- 
hand  side  of  Equation  1  are  real  quantities  obtained  by  measurements; 
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and  the  third  condition  requires  the  continuity  of  the  inverse  function 
Ou.  It  was  felt  for  a  long  time  that  if  at  any  point  u  the  function  Ou 
was  discontinuous,  then  the  solution  f  could  not  be  uniquely  recovered 
from  the  right  hand  side  u.  Hadamard  introduced  the  notion  of  well- 
posedness  by  giving  an  example  of  an  IPP  which  became  a  classical  text¬ 
book  example.  This  example  was  the  famous  Cauchy  problem  for  the  LaPlace 
equation.  Hadamard  did  not  believe  that  an  IPP  represents  any  real 
physical  system.  This  later  conclusion  proved  to  be  erroneous,  and  many 
real  equations  of  mathematical  physics  lead  to  problems  which  are  im¬ 
properly  posed  in  the  sense  of  Hadamard. 

We  now  formulate  an  approach  to  the  question  of  well-posedness  of 
problems  of  the  type  under  consideration.  The  approach  consists  of 
changing  the  notion  of  correctness  by  having  requirements  different  from 
(la),  (lb),  and  (lc)  above.  In  addition  to  the  spaces  U  and  F  and  the 
operator  A,  let  there  be  given  some  closed  set  cpcll .  According  to 
Tikhonov,  the  solution  of  Equation  1  is  properly  posed  if 

(2a)  It  is  "a  priori"  known  that  the  solution  f  exists  for  some 
class  of  data  and  belongs  to  the  given  set  $,  fe$. 

(2b)  The  solution  is  unique  in  a  class  of  functions  belonging  to  $. 

(2c)  Arbitrarily  small  changes  in  u  do  not  carry  the  solution  f 
out  of  $  corresponding  to  arbitrarily  small  changes  in  the  solution  f. 

Upon  denoting  the  image  of  $  after  the  application  to  the  space  F 
of  the  operator  A,  requirement  (2c)-|  can  be  restated  as 

(2c)g  The  solution  of  Equation  1  depends  continuously  on  the  right- 
hand  side  of  u  on  the  set 

If  $  is  a  compact  set  than  according  to  Tikhonov,  if  Equation  1 
satisfies  (2a),  and  (2b),  there  exists  a  function  a  (t),  where  x  is  a 
variable  parameter,  such  that 

(a)  ct(x)  is  a  continuous  nondecreasing  function  with  a(0)  =  0. 


5 


AFML-TR-67-1 21 
PART  VII 

(b)  for  any  fp  fg  e  $  satisfying  the  inequality  p  (Afp  hfg)  £e 

where  p  (\p<p)  is  the  metric  or  measure  of  distance  between  \p  and  <j>  and  e 
is  a  constant,  then  the  following  holds 

p  (fp  i2)<  a(e) 

That  is,  if  a  problem  is  improperly  posed  in  the  metric  spaces  F  and  U, 
it  becomes  properly  posed  in  the  usual  sense  if  F  is  replaced  by  the  sub¬ 
spaces  $  and  $ 

The  reason  for  examining  the  spaces  F,  U  together  with  §,  is  due 
to  the  fact  that  in  real  problems  the  errors  introduced  from  experimental 
measurements  into  the  determination  of  a  set  {u}  usually  result  in  some 
u  being  outside  The  regularization  technique  formulated  by  Tikhonov 
gives  the  possibility  of  constructing  an  approximate  solution  with  a 
certain  guaranteed  degree  of  accuracy  even  though  the  exact  solution  of 
Equation  1  with  approximate  data  either  does  not  exist  or  greatly  deviates 
from  the  "true"  solution. 
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SECTION  III 
QUADRATIC  PROGRAMING 

Consider  the  quadratic  programming  problem  of  finding  {x..}, 
i  =  1,  — ,  n  which  maximizes 


n  n 


b.x.  -  ~  V  V  x.x.g.. 
4-11  24-  4-1  j  ij 

1=1  J=1 


i=l 


subject  to 


X 


c,.x.  1  d 
ki  1  k 


---  ,  m 


(4) 


and  the  non-negativity  conditions 


X.  >  0  i  =  1, - >  n 

1 


(5) 


where  g. .  are  the  elements  of  a  symmetric,  positive  semi -definite  matrix, 

’  J 

i .  e . , 


i] 


=  gii 


(6) 
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and 


n  n 

Z  Z  ekjxjXk  *  0 

j=l  k=l 


(7) 


for  all  x..  It  is  always  possible  to  write  a  quadratic  function  in  the 

J 

form  of  Equation  3  such  that  Equation  6  is  satisfied.  The  restriction 
Equation  7  ensures  that  the  solution  of  Equation  3  is  convex.  There  have 
been  many  algorithms  devised  for  solving  this  problem.  A  few  of  these 
are  one  due  to  Dantzig  (Reference  43),  one  due  to  Thiel  and  Van  de  Panne 
(Reference  44),  another  due  to  Lemke  (Reference  45);  two  based  on  ex¬ 
tensions  of  the  simplex  algorithm  encountered  in  linear  programing  one 
by  Wolfe  (Reference  46)  and  another  by  Beale  (References  47,  48).  In 
addition,  there  are  excellent  review  articles  and/or  books  written  on  the 
details  involved  in  solving  Equations  3,  4,  and  5  (References  49-54). 

In  matrix  notation  Equations  3  through  5  can  be  written  as  maximize, 

B'  X  -  j  X'  GX  (8) 

subject  to 


and 


CX  <  D 


(9) 


X  >  o 


(10) 


where  G  is  positive  semi -definite,  i.e.,  X 1 GX  >_  0  for  all  values  of  X. 
Here  the  "prime"  indicates  the  transpose. 

The  well-known  Kuhn-Tucker  conditions  assert  that  X  is  a  solution  if 
and  only  if  there  exists  a  vector  W  such  that 
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W  >  0 

W'D  -  W'C  X  =  0 


and 


GX  +  C'W  -  B  >  o 


X'GX  +  X'C'W  -  X'B  =  0 

By  making  the  following  substitutions 

V  =  GX  +  C'W  -  B  >  o 


(ID 


(12) 


(13) 

(14) 


(15) 


and 


Y  =  D-  C  X  >  0 


(16) 


the  Kuhn-Tucker  conditions  can  then  be  stated  as  finding  X,  W,  V  and 
Y,  all  >  0,  such  that 


-G  O  E  - 

CEO 


C' 


o 


X 

’-B 

Y 

V 

J 

W  . 

D_ 

(17) 


where  E  are  unit  matrices  and  such  that 


[VW]  ’  y 


V’X  +  W'Y  =  0 


(18) 
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In  the  following  sections  the  method  used  to  incorporate  the  regulari¬ 
zation  technique  of  Tikhonov  into  the  quadratic  programing  scheme  out¬ 
lined  above  will  be  discussed. 
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SECTION  IV 

REGULATION  OF  ILL-POSED  INTEGRAL  EQUATIONS 
OF  THE  FIRST  KIND 

As  an  example  of  the  application  of  Tikhonov's  regularization  tech¬ 
nique,  consider  a  Fredholm  integral  equation  of  the  first  kind, 

Ml 

u(£)  =  Q  [£,f(m)]  =  f  K(f,m)f(m)dm  ,  ^  <  (19) 

M 

o 

Assuming  that  certain  u(£)  functions  exist  which  do  not  have  corresponding 
f(m)  solutions  fulfilling  conditions  (la),  (lb),  and  (1c),  means  Equation 
19  is  an  IPP.  Then  upon  application  of  Tikhonov's  ideas  (Equation  2)  to 
a  special  function  u(?)  there  corresponds  a  solution: 

f(m)  =  R[m,u(£)] 

(20) 

Also  let  an  approximating  function  u(£)  for  u(^)  be  given,  such  that 
||  u  -  u  ||  <6,  where  6  is  known.  It  is  then  required  to  find  f(m),  an 
approximation  to  f(m)  with  an  assigned  precision  ||f  -  f  ||  <_  e  if  6 
is  sufficiently  small.  Letting  MQ  =  0;  =  Mmax;  50=  0  and  ^  =  1; 

assuming  K(£m)  is  continuous  and  if  for  u(£)  =  0  there  exists  just  one 
solution  f(m)  =  0;  then  instead  of  using  the  conventional  functional 
of  calculus  of  variations 


N  [ f(m) ;  u(£)] 


(21) 
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Tikhonov  suggests  application  of  the  functional 


M“  [f(m) ;  u(£)  ]  =  N[f(m);u(f)]  +  a  SI  (n)[f(m)]  (22) 

where  is  the  regularizing  functional 

0(n)(f)  =  f  £  Pi(m)  [f(i)(m)]2|  dm 

4  h=0  1  (23) 

the  P.  (m)  are  positive  continuous  functions,  f^  is  the  i  th  derivative 
with  respect  to  m  and  a  is  an  arbitrary  parameter  which  minimizes  the 
functional  M“. 

Application  of  the  Eulerean  equation  and  applying  boundary  conditions 
results  in 


La  [f] 

n 


a 


dm1 


P.(m) 


£*  i  \ 
dmlJ  j 


M 

max 

K(m,  £  )£(  C  )dC 


(24) 


with  boundary  conditions 


n  (m)  = 


N+I 

£  (-1! 
i=i +  1 


i-i-1 


[P.(m)f(l)(m)](l"i'1) 


M=0,M 


max 


(i  =  1,2,  ,  N  +  1) 
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where 


K(m,  %  ) 


K(£,m)K(£,t,)d£ 


and 


s(m)  = 


K(  £ ,  m)  5(£)d£ 


This  procedure  was  applied  to  a  kernel  of  the  form 


(26) 


(27) 


K(£,s)  =  p  s  e_pS^  /  (l-e"PS)  <28) 


appearing  in  the  theory  of  equilibrium  sedimentation  of  polydisperse 
system,  by  initially  assuming  f(s)  =  const.  Xs  (1-s)  ,  the  set  {u }  was 
computed  using  Equation  19.  m”  of  Equation  22  was  then  minimized  by 
application  of  the  regularizing  technique  resulting  in  the  approximate 
solution  ^(s). Figures  1  and  2  show  the  results  of  the  computation  with¬ 
out  regularization  and  with  regularization,  respectively  (Reference  55). 

During  the  application  of  this  technique  to  a  specific  physical 
problem  it  was  observed  that  when  f(m)  was  multimodal  (bimodal  or  higher) 
then  portions  of  fa(m)  were  negative.  From  physical  considerations  of 
the  problem  of  determining  a  molecular  weight  distribution  from  data 
obtained  from  an  ultracentrifuge  equilibrium  sedimentation  experiment  for 
which  the  kernel  in  Equation  28  is  applicable,  all  fa(m)  should  be 
positive.  Using  these  considerations  regularization  was  incorporated 
into  (LP)  linear  programing  (Reference  56)  using  Dantzig's  Simplex 
algorithm  (Reference  57).  The  regularized  LP  technique  gave  good  results 
up  through  a  trimodal  distribution.  For  higher  multimodal  distributions 


13 


AFML-TR-67-1 21 
PART  VII 


the  computed  fa(m)  were  very  erratic  and  the  computational  error  was 
large.  However,  since  the  functional  to  be  minimized  (Equation  22)  is 
quadratic,  it  seemed  only  natural  to  apply  quadratic  programing. 

In  the  following  the  incorporation  of  regularization  into  the  quadratic 
programing  algorithm  given  by  Boot  (Reference  54)  is  discussed. 
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SECTION  V 

INCORPORATION  OF  REGULARIZATION 
INTO  QUADRATIC  PROGRAMING 

In  all  applications  involving  the  kernel  given  by  Equation  28  it  has 
been  found  that  sufficiently  satisfactory  results  were  obtained  when  the 
P.j (m)‘s  in  Equation  23  were  equated  to  constants.  Therefore  the  functional 
in  Equation  22  to  be  minimized  was  restricted  to  become 


Ma  [f(m) 
n 


uU)]  =  N  [f(m)  ;  u(£)]  +X!an 

n 


f!,n)  [f] 


(29) 


where  now 


&(n)[f] 


M 

/max  ,  .  , 

[f  n  (m)]  dm 


(30) 


f^n^(m)  being  the  n^  derivative  of  f(m)  with  respect  to  m^(f^(m)  = 
dnf(m)/dmn).  The  n*^  derivative  of  f(m)  for  n  =  1,  2,  3,  .  .  .  can  be 
approximated  by  various  numerical  techniques.  In  this  specific  case, 
assume  h  is  the  constant  increment  associated  with  the  mesh  for  m.  Then 
f|.n^  can  be  approximated  by 

J 


£(n)  =  1  £  (TVi-i)"  f. 

3  hn  k^O  J“P+k 


(31) 


where  ^  are  the  binomial  coefficients;  and  p  =  n  for  n  odd;  and  p  =  n-1 
for  n  even.  Then  Equation  30  becomes 


“"['I'Fin!  |  («)(;)(-') 

n  j  =  i  i  =  i  k=o  £-o 


£*k 


fj-P+kVp+^  (32) 
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which  in  matrix  notation  will  be 

r  ■)  ,  .(n) 

Q.  [f J  *  f1  A  f 


(33) 


where  A  is  a  matrix  whose  elements  are  zero  except  for  diagonal  and  near 
off  diagonal  elements  for  which  if  r  =  j  -  p  +  k  and  s  =  i  -  p  +  1  (as  in 
Equation  32)  then 


(34) 


with  the  boundary  conditions  1  £  r  £  J  and  for  s  <  1 ,  then  s  =  |s|  +  1  or 
s  >  J,  then  s  =  2J  -  s  +  1 . 

Next  consider  Equation  21.  Let  us  express  this  in  matrix  notation, 
where  as  in  Equation  1  the  operator  (kernel  multiplied  by  appropriate 
integration  constants  for  numerical  evaluation)  will  be  designated  by  A 
=  {a..};  u  =  (u . >  i  =  1,2,...,I;  and  f  =  {f .}  j  =  l,2,...,d.  Thus 

I J  1  J 

Equation  21  can  be  expressed  as 


fWAf  -2u'Af +  u'u 


(35) 


Neglecting  the  last  term  in  Equation  35  which  is  a  constant,  and  using 
the  result  of  Equation  33,  the  functional  of  Equation  29  expressed  in 
matrix  notation  is 


a 

Mn 


f'A'Af  -  2u'Af  +  £anf 


f 


n 


(36) 
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or  upon  dividing  by  2  to  be  in  correspondence  with  Equations  3  and  8  and 
multiplying  by  (-1),  the  functional  to  be  maximized  will  be 


=  u'Af  -  4-f1  [  A'A+Xan  A1"*  1  f 

c  L  n  J  (37) 


J 

(Equation  8),  subject  to  £  t.  f.  <_  const,  and  all  f.>  0. 

J  J  J 

j=l 

This  is  now  a  suitable  quadratic  program  formulation.  In  the 
following  a  particular  kernel  will  be  used  and  a  computer  simulation 
experiment  where-in  analogous  experimental  data  {u  :  u  =  Af}  is  generated 
from  an  assumed  set  {7}  and  the  back  solution,  determining  f  from  U  will 
be  discussed.  Since  in  a  real  experimental  situation  the  original  f 
would  be  unknown,  f  and  f  are  presented  only  for  illustrative  purposes. 
All  computations  were  performed  so  as  to  choose  that  set  of  an’s  (usually 
a  single  an  sufficed)  which  yielded  a  minimum  for  | |U  -  u||  i.  e.,  the 

error  criterion  was  to  choose  that  {f}  in  correspondence  with  inf  { | | U 

0(.n 
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SECTION  VI 

APPLICATION  AND  RESULTS 


The  first  step  in  proving  the  utility  of  Equation  37  was  to  establish 
a  kernel  which  represented  an  IPP  in  a  real  physical  situation.  Such  an 
expression  is  Equation  28.  The  computational  work  was  then  related  to 
the  following  integral  equation  of  the  first-kind. 


u 


f  (m)  dm 


(38) 


-5 

where  8  =  const,  and  0  <_  E  <_  1 .  In  all  cases  8  =  4  x  10  .  For  ummodal 

and  trimodal  distributions,  f(m),  a  41-point  mesh  was  used  for  E  and  m; 
for  a  tetramodal  distribution  a  51-point  mesh  and  for  an  asymmetrical 
bimodal  and  pentamodal  a  59-point  mesh  was  used.  That  is,  if  N  equals 

n  j 

the  number  of  intervals  in  our  mesh  then  £  =  TN  where,  n-|  =  1,  2,  ..., 

N  -  1  and  m  =  n2mniax/N,  n^  =  1 ,  2,  . . .  ,N  -  1 .  All  integrations  were 
performed  using  Simpson's  quadrature  formula  for  equidistant  points.  It 
was  felt  that  in  real  problems  this  would  be  sufficient  and  it  was  not 
the  purpose  of  this  research  to  study  how  to  minimize  machine  round-off 
errors. 

An  initial  functional  distribution  f(m),  unimodal  through  pentamodal, 
was  assumed.  Then  Equation  38  was  used  to  compute  a  set  of  values  for 
lT(E).  These  were  then  assumed  to  be  our  experimental  values. 


Next,  quadratic  programing  with  regularization  was  applied  (Equation 

■  ^  ~QL 

37).  For  a  given  a  ,  the  corresponding  set  {f}  which  minimized  was 

n  c Xfl  n 

computed.  Then  through  application  of  Equation  38  the  corresponding  set 
(u  (£)}  was  evaluated.  The  a  which  yielded  inf  {Mu"  -  u 1 1 }  was  the 
final  caused.  Further  searching  for  an  c*n  with  more  significant  digits 
would  have  decreased  the  error  analysis  criterion  but  for  our  purposes 
two  significant  figures  were  considered  satisfactory.  Finally  the 
initial  f(m)  and  the  f^m)  were  plotted  in  order  to  compare  the  distri¬ 
butions. 
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These  distributions  along  with  the  computed  data  are  presented  in 
Figures  2  through  7  and  Tables  I  through  V,  respectively.  To  show  the 
need  for  regularization  some  figures  are  presented  with  the  results 
obtained  when  no  regularization  -  only  quadratic  programing  was  used.  In 
addition  it  should  be  noted  that  the  fewer  points  per  mode  the  less  the 
precision.  This  is  especially  noticeable  when  comparing  the  unimodal  and 
pentamodal  distributions.  In  the  former,  41  points  were  used  per  mode 
whereas  in  the  latter  there  were  only  about  11  or  12  points  per  mode. 

Due  to  round-off  errors,  storage  space  in  a  high  speed  digital  computer, 
and  computational  time  the  present  computation  was  limited  to  using  no 
more  than  about  11  points  per  mode  for  the  pentamodal  distribution.  As  a 
demonstration  of  this  necessity  to  sample  a  sufficient  number  of  molecular 
weights,  the  following  test  was  performed.  Starting  with  nine  molecular 
weights  the  initial  f  (m*)  was  computed.  From  these  functional  values  the 
corresponding  set  u  (£*)  was  inferred  on  a  41-point  mesh.  This  number 
of  values  was  used  to  compute  the  corresponding  set  fa(m*),  in  the  same 
fashion  as  f(m*).  Finally,  the  set  fa(m*)  was  used  to  compute  an  anal¬ 
ogous  set  u  (£*).  Figure  8  shows  a  comparison  of  u"  (£)  computed  from  the 
f  (m)  with  41  molecular  weight  (Figure  3)  with  IT  (£*)  and  u  (£*)  computed 
using  a  nine-point  molecular  weight  mesh. 
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SECTION  VII 
CONCLUSIONS 

The  need  for  knowing  the  molecular  weight  distribution  of  synthetic 
polymers  first  led  the  authors  to  the  ill-posed  inverse  problem  associated 
with  Equation  38.  Scientists  have  investigated  the  feasibility  of  this 
determination  for  the  past  30  to  40  years.  All  types  of  well  founded 
mathematical  theories  were  applied,  but  each  would,  in  general,  only  be 
applicable  for  specific  types  of  distributions.  It  was  only  recently 
realized  that,  instead  of  apologizing  for  the  kernel  of  Equation  38  being 
ill-conditioned,  the  entire  problem  was  mathematically  ill-posed  in  the 
Hadamard  sense.  It  should  be  challenging  to  derive  another  expression 
for  determining  a  molecular  weight  distribution  from  equilibrium  sedi¬ 
mentation  data  which  might  be  a  well-posed  problem.  Meantime,  (since 
time  and  economics  prevented  such  a  diversion)  application  of  Tikhonov's 
technique  of  regularization  has  enabled  reliable  results  to  be  obtained. 
Good  results  were  obtained  for  unimodal  through  tetramodal  distributions. 
Poor  results  were  obtained  for  a  pentamodal  distribution.  The  results 
indicate  that  even  if  the  experimental  data  u(c)  are  precise,  a  "poor 
fit"  MWD  will  be  obtained  if  the  sampling  size  of  molecular  weights  is 
too  small.  It  can  be  estimated  that  a  lower  limit  on  the  number  of 
molecular  weights  per  mode  or  per  peak  has  to  be  about  20  in  order  to 
obtain  a  good  "fit".  Ten  molecular  weights  per  peak  gave  poor  results. 

To  assure  such  a  good  "fit"  a  bimodal  distribution  would  require  a  40- 
point  mesh  minimum.  Unfortunately,  because  of  the  computer  storage 
limitations,  as  well  as  an  extensive  computation  time,  the  mesh  could 
not  exceed  61  points.  This  number  was  adequate  to  compute  a  trimodal , 
barely  adequate  to  compute  a  tetramodal  distribution,  and  entirely  in¬ 
adequate  to  compute  a  pentamodal  distribution.  Considering  these 
limitations,  the  computation  of  higher  multimodal  distributions  were  not 
attempted. 

In  addition,  a  larger  molecular  weight  mesh  would  also  require  a 
corresponding  larger  number  of  discernible  u  (?).  For  the  ultracentri¬ 
fugal  techniques  this  would  require  the  use  of  longer  column  lengths  for 
solutions  investigated. 
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APPENDIX 

QUADRATIC  PROGRAM  ALGORITHM 


Computer  programs  were  written  to  solve  the  following  problems* Find 
the  values  of  x-j ,  x2>*-’>xn  that  maximize 


subject  to 


where 


11 

°  1 2 

c 

In 

c 

c 

21 

22 

2n 

c 

c 

kl 

k2 

kn 

x.  >  O  i  =  1,2, 


b 

b 

b  1 

11 

12 

In 

b 

b 

.  .  .  b 

21 

22 

2n 

b  , 

b 

b 

nl 

n2 

’  ’  nn 

is  positive  semi-definite 


In  matrix  form  we  have: 

Find  the  value  of  *  that  maximizes 

A'x  -  1/2  x'  m 
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subject  to 

CX<D 
X>0 

where  B  is  positive  semi -definite,  i.  e.,  X'BX>_0  for  all  values  of  X, 
This  problem  can  be  reformulated  by  introducing  k  non-negative  slack 
variables  (y-j 
Find  the  values  of  X,  Y  that  maximize 

[?]  -  [l  §]  K 

subject  to 

[C  i]  [*] .  D 
X  >  o 
Y  >  ° 

Using  the  Kuhn-Tucker  conditions,  it  can  be  shown  that  [X  Y]/  where 
prime  denotes  transpose,  is  the  solution  to  this  problem  if  and  only  if 

1.  [X  Y]'  >  0 

2.  There  exists  a  vector  [V  W]  of  non-negative  elements  such  that 

[v  w]7  j^J  =  V'X  +  W’Y  =  o 

3.  The  vectors  [X  Y];  and  [V  W]/  satisfy  the  system  of  linear  equations 


y^)  =  Y  (Reference  53),  and  stating  the  problem  as 


Dantzig's  alogrithm  as  presented  by  Boot  (Reference  53)  is  used.  Th 
procedure  begins  with  a  basic  feasible  standard  form  solution  (X  Y  V  W)‘ 
(0  D  -AO)'  of  the  system  of  linear  equations  above.  The  system  of 
linear  equations  has  m  +  k  equations  and  2(m  +  k)  unknowns. 
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A  basic  solution  is  a  solution  determined  by  setting  m  +  k  of  the  var¬ 
iables  equal  to  zero  and  solving  the  remaining  variables.  A  basic  feasible 
solution  is  a  basic  solution  that  has  only  non-negative  values  for  [X  Y]. 
Standard  and  nonstandard  basic  feasible  solutions  are  defined  as  follows. 
Let  Z1  =  [X  Y] 1  and  U1  =  [V  W]1.  If  a  basic  feasible  solution  is  such 
that  no  pair  of  corresponding  Z  and  U  variables  consist  of  two  nonzero 
elements,  the  solution  is  in  standard  form,  otherwise  the  basic  feasible 
solution  is  in  nonstandard  form. 

In  a  basic  feasible  solution  of  the  system  of  linear  equations,  the 
variables  that  are  set  equal  to  zero  are  called  nonbasis  variables,  the 
remaining  are  called  basis  variables  and  comprise  the  basis.  The  algo¬ 
rithm  consists  of  adding  a  variable  to  the  basis  and  deleting  a  variable 
from  the  basis.  This  is  better  explained  by  writing  the  system  of  linear 
equations  as  a  linear  combination  of  vectors  equal  to  a  vector. 

Let  P  equal  the  m™  column  of  the  matrix 
m 


-B 

O 

I 

-C 

C 

I 

o 

o 

and  let  P*  =  [-A  D],  Then  the  system  of  linear  equations  can  be  written 

Z1  P1  +  Z2P2  +  +  Zn+k  +  UlPn+k+l  +  U2Pn+k+2  +  \**  +  Un+kP2N+2k  =  Po 

Let  P  ,  m  =  1 ,2. . . ,  n  +  k  be  the  values  of  the  basis  variables  and  let 
0m  =  the  subscript  of  the  associated  P  vector  for  the  m  th  basis  variable 
m  =  1 ,2, . . .  ,n  +  k 

The  rules  for  adding  a  variable  to  the  basis  are: 

1.  If  the  basic  feasible  solution  is  in  standard  form,  that  particular 
non-basic  Z-variable  should  enter  the  basis  whose  corresponding  Uh  has 
(in  absolute  value)  the  largest  negative  P^. 

2.  If  the  basis  feasible  solution  is  nonstandard  and  (Z^»  U^)  is  the 
nonbasic  pair,  then  should  enter  the  basis. 
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Let  the  P  vector  corresponding 
basis  be  represented  by  (T.  L,  ... 
variable  from  the  basis  are: 


to  the  variable  that  is  to  enter  the 
Tn+|<)1-  The  rules  for  deleting  a 


1.  If  the  basic  feasible  solution  is  standard,  let  be  the  variable 

that  is  to  enter  the  basis.  Find  the  value  of  m  that  corresponds  to  the 

smallest  positive  ratio  P  /T  while  only  considering  those  m's  such  that 
r  m  m  J  3 

jme  {1 ,2,...,n  +  k,  n  +  k  +  h}. 


2,  If  the  basic  feasible  solution  is  nonstandard,  let  (Z^,  U^)  corre¬ 
spond  to  the  pair  that  are  both  basic.  Find  the  value  of  m  that  corre¬ 
sponds  to  the  smallest  positive  ratio  pm/Tm  while  only  considering  those 
m's  such  that  jme  {1 ,2,...,n  +  k,  n  +  k  +  h} 


The  algorithmic  recycling  is  terminated  when  all  of  the  basic 
variables  are  nonnegative,  i.  e.,  when  Pm  >  0;  m  =  l,2,...,n  +  k. 
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TABLE  I 


COMPUTATIONAL  RESULTS  FOR  THE  UNIMODAL  DISTRIBUTION 


FROM  A  41-POINT  MESH 


No. 

m 

f(m)»  106 

?"(m)x  106 

u{f  ) 

u(f  ) 

1 

3,571 

0.108 

0.200 

3.2052 

3.2061 

3 

10,714 

0.880 

0.791 

2.7151 

2.7159 

5 

17,857 

2.200 

2. 170 

2.3069 

2. 3075 

7 

25,000 

3.858 

3.931 

1.9660 

1 . 9664 

9 

32,143 

5.670 

5.737 

1.6806 

1.6809 

11 

39,286 

7.474 

7.460 

1.4410 

1.4413 

13 

46,429 

9.135 

9.061 

1.2393 

1.2396 

15 

53,  571 

10.542 

10.477 

1.0692 

1.0694 

17 

60,714 

11.609 

11.623 

0.9252 

0.9254 

19 

67,857 

12.274 

12. 356 

0.8031 

0.8033 

21 

75,  000 

12.500 

12.578 

0.6992 

0.6993 

23 

82, 143 

12.274 

12.276 

0.6106 

0.6107 

25 

89,286 

11.609 

11.523 

0.5348 

0.5349 

27 

96,429 

10.542 

10.438 

0.4698 

0.4699 

29 

103,  570 

9. 135 

9.095 

0.4139 

0.4140 

31 

110,710 

7.474 

7.554 

0. 3657 

0.3658 

33 

117,860 

5.670 

5.813 

0.  3240 

0.  3241 

35 

125,000 

3.858 

3.912 

0.2879 

0.2880 

37 

132, 140 

2.200 

2.076 

0.2565 

0.2566 

39 

139,290 

0.880 

0.761 

0.2291 

0.2292 

41 

146,430 

0.108 

0.270 

0.2052 

0.2053 
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TABLE  II 

COMPUTATIONAL  RESULTS  FOR  THE  ASYMMETRICAL  BIMODAL 
DISTRIBUTION  FROM  A  59-POINT  MESH 


No. 

m 

f(m)  x  10^ 

Ta(m)  x  106 

uif  ) 

u(£  ) 

1 

2,  500 

0.  173 

o.  no 

2.8642 

2.8644 

3 

7,  500 

1.388 

1.  565 

2.5625 

2.5727 

5 

12,  500 

3.401 

3.639 

2.3163 

2.3164 

7 

17,500 

5.834 

5.484 

2.0906 

2.0908 

9 

22, 500 

8.259 

7.  587 

1.8915 

1.8916 

11 

27, 500 

10.706 

10.483 

1.7154 

1.7156 

13 

32, 500 

12.656 

13. 181 

1.5594 

1.5595 

15 

37,500 

14.047 

15.050 

1.4208 

1.4209 

17 

42,500 

14.769 

15.621 

1.2975 

1.2976 

19 

47,500 

14.769 

14.808 

1.1875 

1.1876 

21 

52, 500 

14.047 

13.408 

1.0891 

1. 0892 

23 

57,500 

12.656 

11.481 

1.0010 

1.0010 

25 

62,500 

10.706 

9.538 

0.9218 

0.9216 

27 

67,  500 

8.359 

7.974 

0.8506 

0.8506 

29 

72, 500 

5.834 

6.434 

0.7863 

0.7863 

31 

77,  500 

3.520 

4.863 

0.7282 

0.  7282 

33 

82, 500 

2.  316 

3.853 

0.6756 

0.6756 

35 

87, 500 

2.385 

2.926 

0.6278 

0.6278 

37 

92, 500 

3.670 

2.  530 

0.  5844 

0.  5844 

39 

97,500 

5.057 

3.567 

0.  5448 

0.  5448 

41 

102,500 

6. 184 

5.749 

0.  5086 

0.  5086 

43 

107, 500 

6.914 

7.600 

0.4755 

0.4755 

45 

112, 500 

7. 167 

7.632 

0.4452 

0.4452 

47 

117,500 

6.914 

6.817 

0.4173 

0.4173 

49 

122,  500 

6. 184 

6. 124 

0. 3917 

0. 3917 

51 

127, 500 

5.057 

4.902 

0.3681 

0.  3681 

53 

132, 500 

3.670 

3.765 

0.  3463 

0.346  3 

55 

137,  500 

2.212 

2.483 

0. 3261 

0.3261 

57 

142,  500 

0.929 

0.964 

0.  3075 

0.  3074 

59 

147,  500 

0. 119 

0.000 

0.2902 

0.2901 

=  7.9  x  10'7 

II  u  -  vT  II  =  8.73x10^ 

=  5.6  x  10"10 

£  =  Number /6l 
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TABLE  IV 

COMPUTATIONAL  RESULTS  FOR  THE  TETRAMODAL  DISTRIBUTION 
FROM  A  53-POINT  MESH 


No. 

m 

f(m)  x  10^ 

7\m)  x  106 

u  U  ) 

«(£) 

1 

2,778 

0.592 

0.  000 

3.2556 

3.2556 

3 

8,  333 

3.991 

3.  374 

2.8297 

2.8297 

5 

13,889 

7.910 

11.405 

2.4678 

2.4678 

7 

19,444 

10.  329 

12.895 

2. 1595 

2. 1595 

9 

25,  000 

10.251 

6.929 

1.8964 

1.8963 

11 

30,  556 

7.706 

0.000 

1.6712 

1.6711 

13 

36, 111 

3.807 

2.223 

1.4780 

1.4780 

15 

41,667 

3.064 

8.468 

1.3120 

1.3119 

17 

47,222 

6.616 

13.278 

1.1688 

1. 1688 

19 

52,778 

9.731 

13.488 

1.0450 

1.0450 

21 

58,333 

10.583 

9.582 

0.9378 

0.9378 

23 

63,889 

8.841 

3.  356 

0.8446 

0.8445 

25 

69,444 

5.199 

0.  000 

0.7633 

0.7633 

27 

75,000 

2.754 

0.593 

0.6923 

0.6923 

29 

80,556 

5.199 

3.  511 

0.6300 

0.6300 

31 

86, 111 

8.841 

11.522 

0.  5752 

0.5753 

33 

91,667 

10.583 

16.024 

0.5269 

0.  5270 

35 

97,222 

9.731 

14.384 

0.4842 

0.4842 

37 

102,780 

6.616 

7.620 

0.4463 

0.4463 

39 

108, 330 

3.064 

1.743 

0.4126 

0.4126 

41 

113,890 

3.807 

0.000 

0. 3824 

0.3825 

43 

119,440 

7.706 

1.256 

0.3555 

0.3555 

45 

125,  000 

10.251 

9.745 

0. 3312 

0.3313 

47 

130, 560 

10.  329 

16.270 

0. 3094 

0.3095 

49 

136, 110 

7.910 

10. 028 

0.2897 

0.2898 

51 

141,670 

3.991 

2.  310 

0.2718 

0.2719 

53 

147,220 

0.  592 

0.000 

*2.556 

0.2557 

q2 

2. 1  x  10‘7 

II  u  -  u 

||  r  4.  0120  X 

io-5 

a  =  7.0x10  |  =  Number/55 
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TABLE  V 


COMPUTATIONAL  RESULTS  FOR  THE  PENTAMODAL  DISTRIBUTION 
FROM  A  59-POINT  MESH 


No. 

m 

f  (m)  x  10^ 

7a(m)  x  106 

u(  f ) 

1 

3,  000 

0.593 

2.609 

3.8132 

3.8089 

3 

9,000 

3.853 

2.747 

3.2709 

3.2683 

5 

15,000 

7.245 

8.257 

2.8166 

2.8153 

7 

21,000 

8.742 

12.430 

2.4351 

2,4347 

9 

27,000 

7.606 

9.559 

2. 1138 

2.  1141 

11 

33,000 

4.391 

5.370 

1.8425 

1.8433 

13 

39,000 

2.281 

2.880 

1.6128 

1.6139 

15 

45,  000 

4.922 

3.231 

1.4177 

1.4190 

17 

51,000 

7.923 

7.604 

1.2516 

1.2529 

19 

57,000 

8.687 

4.621 

1.1096 

1. 1110 

21 

63,000 

6.843 

4.082 

0.9880 

0.9892 

23 

69,000 

3.329 

6.802 

0.8834 

0.8845 

25 

75,000 

2.  592 

10.928 

0.7931 

0.7941 

27 

81,000 

5.935 

5.983 

0.7150 

0.7159 

29 

87,000 

8.410 

0.000 

0.6471 

0.6478 

31 

93,000 

8.410 

4.985 

0.  5879 

0.5885 

33 

99,000 

5.935 

11.253 

0.5362 

0.5366 

35 

105,000 

2.592 

20.421 

0.4907 

0.4910 

37 

111,000 

3.329 

1.102 

0.4507 

0.4509 

39 

117,  000 

6.843 

0.000 

0.4153 

0.4154 

41 

123,000 

8.687 

0.000 

0.3839 

0.3838 

43 

129,000 

7.923 

0.000 

0.  3560 

0. 3558 

45 

135,000 

4.922 

14.118 

0.3310 

0.3308 

47 

141,000 

2.281 

4.616 

0.  3087 

0.3083 

49 

147, 000 

4.391 

1.905 

0.2886 

0.2881 

51 

153,  000 

7.606 

14.409 

0.2704 

0.2699 

53 

159,000 

8.742 

7.154 

0.2540 

0.2534 

55 

165,000 

7.245 

0.000 

0.2391 

0.2385 

57 

171,000 

3.853 

0.000 

0.2256 

0.2249 

59 

177,000 

0.593 

7.073 

0.2132 

0,2125 

34 


Figure  2.  Unimodal  Distribution  by  Variational  Calculus  with 
Regularization 
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Figure  5.  Symmetrical  Trimodal  Distribution,  Solid  Line  the  Original 
Distribution,  the  Circles  Quadratic  Programing  and  Regular¬ 
ization,  the  Histogram  to  1/10  Scale  Using  Only  Quadratic 
Programing 


Histogram  is  zD  scale 
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Figure  7.  Symmetrical  Pentamodal  Distribution,  Solid  Line  the 

Original  Distribution,  the  Circles  Quadratic  Programing 
with  Regularization,  the  Histogram  to  1/10  Scale  Only 
Quadratic  Programing 


if  the  Computed  Curve  u  (€*)  and  IT  (£*)  from 
r  (C)  as  a  function  of  S 


PROGRAM  LISTING 
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******  ******  ****** 


PROGRAM  REQUAD 
PURPOSE 

PROGRAM  REAOS  ALPHA (NLAST)  AS  DATA  FROM  NFIRST 
THROUGH  NCOOE.  THEN  PROGRAM  CONTINUES  FROM 
NCODEtl  THROUGH  NUPP  SEARCHING  FOR  MINIMUM  FOR 
EACH  DERATIVE  RETAINING  PREVIOUS  VALUES.  IF 
NCODE  =  0  ,  SEARCH  BEGINS  WITH  NFIRST, 

IF  NFL AG. GT • 0  PROGRAM  READS  ONE  VALUE  OF  ALPHA  AND 
COMPUTES  FOR  ONLY  THIS  ONE  VALUE 


USAGE 

PROGRAM  REQUAD (TAPE5, OUTPUT, TAPE6*0UTPUT) 


PROGRAM  REQUAD (TAPES, OUTPUT, TAPE6=0UTPUT) 

COMMON/ZYT/U, XK,S,X,Z,BK, B,ZP, A, ALPHA, R, XX, UP, OLX,DLS, 
*IMAX,NMAX,NF 

1IRST, NLAST, FACTO, TCOST,TPIV 
C  DIMENSIONS  FOR  COMMON 

DIMENSION  U(60) ,XK(60,6Q) ,S<60) ,X(60) ,Z(60) ,BK (60,601 , 
*B(60),ZP(60) 

1, A (60, 60) , ALPHA (10) ,R(60> ,XX(60) ,UP(60) 

DIMENSION  I3ASIS(60) , RESULT <601 
C 

READ (S, 200) 

MRITE(6,1Q00) 

MRITE(6»200) 

MRITE(6,100i) 

C  NMAX  AND  IMAX  MUST  BE  ODD  INTEGERS 

READ  (5,201)  NNM AX, UMAX, Ml, M2 
NMAX  =  NNMAX  -  2 
LMAX  =  NMAX  ♦  1 
IMAX  =  UMAX 

DLS  =  FLOAT (NNM AX-1) / ( H2-M1) 

DLX  =  FLOAT (IIMAX-1) 

DO  1  I  =  1 , IMAX 
X(I)  =  FLOAT ( I-l) /DLX 
1  CONTINUE 
COF  =  0. 

DO  2  I  =  1, NMAX 

S(I)  =  Ml  t  FLOAT ( I ) /DLS 

FI  =  0, 

IF(S(I).GT.W2.0R.S(I),LT.M1)  GO  TO  53 
A1  =  (S(I)  -  Ml)**2 
A2  =  (S(I)  -  M2)**2 
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Fi  =  A1*A2 
53  ZCI)  =  Fi 
KNUM  =  1/2 
JNUM  =  <I+l)/2 
IF(JNUM.NE.KNUM)  60  TO  5i 
SI6  =  2. 

GO  TO  52 

51  SIG  =  4. 

52  COF  =  COF  +  SIG*Z(I)/(3.*DLS) 

2  CONTINUE 

DO  5  I  =  i,NMAX 
5  Z(I)  =  ZID/COF 

G  XSIG  =  LAMBDA  IN  THE  THEORY,  SEE  FUJITA'S  EQUATION 

READ(5, 103)  XSIG 
00  4  I  =  1, IMAX 
COEF  =  0. 

DO  3  J  *  1, NMAX 
Ai  =  XSIG*S(J) 

A2  =  Al 

A3  =  EXP(-Ai*X(I)) 

A4  =  EXP(-Al) 

A5  =  1.  -  A4 

A6  =  A2*A3/A5 

XKCI , J)  =  A5 

KNUM  =  J/2 

JNUM  =  (J+D/2 

IF(JNUM.NE.KNUM)  GO  TO  42 

SIG  =  2. 

GO  TO  43 

42  SIG  =  4. 

43  COEF  =  COEF  +  SIG*XK< I , J) *Z< J) / <3.*DLS> 

3  CONTINUE 

C  CALCULATION  OF  UCZI)  BY  SIMPSON'S  FORMULA 

4  UCI)  =  COEF 
CALL  REG2 

81  CONTINUE 

READ <5, 101)  NCODE,NFIRST, NUPP, NFL AG, FACTO, TCOST,TPIV 

DO  30  NLAST  s  NFIRST, MUPP 

IF(NLAST.LE.NCODE)  GO  TO  31 

IF(NFLAG.GT.O)  GO  TO  40 

READ (5,100)  LPHA1,LPHA2 

LXP  =  IABSILPHA2  -  LPHA1)  ♦  1 

NUM  =  0 

00  20  II  =  1,  LXP 
IXP  =  LPHAi  +  II  -  1 
DO  21  KL  =  1,9 

ALPHA (NLAST)  =  FLOAT(KL)*10.**IXP 
C  OBTAIN  MODIFIED  MATRIX 

CALL  REG3 
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C  OBTAIN  INVERSE  SOLUTION 

CALL  QUADi (RESULT* IBASIS) 

DO  33  I  =  l.LMAX 
33  ZP(I)  =  0. 

DO  32  I  *  i,LMAX 
J  =  IBASIS ( I ) 

IF(J.GT.LMAX)  GO  TO  32 
ZP(J)  =  RESULT (I) 

32  CONTINUE 

C  EVALUATE  ERROR 

CALL  REGMUAVG) 

WRITE (6,4000)  ALPHA(NLAST) ,UAVG 
4000  FORMAT (1H  ,iP2E12.5) 

IF(KL.EQ.i.AND.II.EQ.l)  GO  TO  22 
IFCUAVG.GE.AVGil  GO  TO  28 
AVGi  =  UAVG 

C  STORE  MINIMUM  ERROR  AND  CORRESPONDING  ALPHA 

XM  =  ALPHA (NL AST) 

NX  =  IXP 

NUM  =  0 
GO  TO  21 

22  AVGI  *  UAVG 

S  STORE  FIRST  ALPHA  USEO  AND  ASSOCIATED  ERROR 

XM  =  ALPHA(NLAST) 

NX  =  IXP 

GO  TO  21 

28  CONTINUE 
21  CONTINUE 
20  CONTINUE 

23  CONTINUE 

IF(XH.EQ.10.**NX)  GO  TO  60 
XMM  *  XM  -  10.**NX 
GO  TO  62 

60  XMM  =  9.*10.**(NX-1) 

62  CONTINUE 

DO  25  I  =  1,20 

ALPHA(NLAST)  =  XMM  +  FLOAT (1-1) *10 .**( NX-i) 
CALL  REG 3 

CALL  QUADI (RESULT, IBASIS) 

DO  55  K  =  1 , LMAX 

55  ZPOO  =  0. 

DO  56  K  =  1, LMAX 
J  =  IBASIS ( K) 

IF(J.GT.LMAX)  GO  TO  55 
ZP(J)  =  RESULT (K> 

56  CONTINUE 

CALL  REGMUAVG) 

IF(I.EQ.i)  GO  TO  26 
IF (UAVG* GE» AVGI)  GO  TO  25 
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AVG1  =  UAVG 

XM  =  ALPHA(NLAST) 

GO  TO  25 

26  AVGi  =  UAVG 

XM  =  ALPHA(NLAST) 

25  CONTINUE 

27  ALPHA (NLAST)  =  XM 
GO  TO  61 

40  REAO (5, 102)  ALPHA (NLAST) 

IF  COMPUTATION  PROCEEDS  FOR  ONLY  ONE  ALPHA  BEGIN 
HERE 

61  CONTINUE 

OBTAIN  MODIFIED  MATRIX 
CALL  REG3 

OBTAIN  INVERSE  SOLUTION 
CALL  QUADi (RESULT *  I BA 5 IS) 

DO  34  I  =  1 » LMAX 

34  ZP(I)  =  0. 

DO  35  I  =  1,LMAX 
J  =  IBASIS ( I) 

IF(J.GT.LMAX)  GO  TO  35 
ZP(J)  =  RESULT (11 

35  CONTINUE 

EVALUATE  ERROR 
CALL  REG4 (UAVG) 

DO  24  I  =  1 ,  IMAX 

Z  <11  =  ORIGINAL  DISTRIBUTION 
ZP(I)  =  BACK  SOLUTION 

U(I>  =  CORRESPONDS  TO  INPUT  DATA,  COMPUTED  USING 
Z(I) 

UP(I)  =  BACK  SOLUTION  COMPUTATION 

S(I)  =  VARIABLE  FOR  Z(I>,  CORRESPONDING  TO 

MOLECULAR  HEIGHT 

WRITE (6, 20 01)  I,ZP(I) , I,Z(I) ,I,UP(I) ,I,U(I) ,I,S(I) 
24  CONTINUE 

WRITE(6,2000)  <  I ,  ALPH  A  ( I)  ,  I  =  NFIRST, NLAST) 
HRITE(6,2002)  UAVG 
WRITE(6,104)  XSIG 
WRITE(6,1000) 

GO  TO  30 

31  REAO (5,102)  ALPHA (NLAST) 

30  CONTINUE 

READ (5,100)  I T1 , IT2 
IF(ITl.EQ.O)  GO  TO  99 
GO  TO  81 
99  CONTINUE 

WRITE(6»1001) 

WRITE(6,7000) 

STOP 
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100  F0RMAT(2I3) 

101  FORMAT  (4I2,1P3E8»1) 

102  FORMAT (1PE14*7) 

103  FORMAT (E10 • 3) 

10%  FORMAT (1H  ,7HXSIG  =  ,1PE10.3) 

200  FORMAT (72H 

1  ) 

201  F0RMAT(2I2,1P2E12.5) 

300  FORMAT (10F8*  5) 

1000  FORMAT (1H1) 

1001  FORMAT (//) 

2000  FORMAT (1H  ,SHALPHA(,I2,4H)  =  ,1PE14.7) 

2001  FORMAT (1H  ,7HZ-CALC(,I2,4H)  =  ,E12.5,2X,7HZ-TRUE( ,12,4 
*H)  =  ,E12.5, 

12X,7HCALC  U(,I2,4H)  =  ,E1 2. 5, 2X, 2HU( , 1 2, 4H)  =  ,E12.5,2 

*X,2HS(,I2,4H 

2)  =  ,E12.5) 

2002  FORMAT (1H  ,30HSQRT  OF  SUM  (UP ( I ) -U (I) ) **2  *  ,E12.5) 
7000  FORMAT (1H  , 20X, 6C5X,10HEND  OF  RUNI71H1) 

END 

*MM****»M**M*f*»»M»*»*******»M**M,*«**y,*,Mm 


SUBROUTINE  REG4 (UAVG) 

PURPOSE 

THIS  SUBROUTINE  PROCESSES  THE  COMPUTED  ZP(I), 
CALCULATES  UP(I>  AND  THE  ERROR  CRITERION 

USAGE 

CALL  REG4CUAVG) 


SUBROUTINE  REG4 (UAVG) 

COMMON/ZTT/U, XK,S,X,Z,BK,B,ZP , A, ALPHA ,R,XX,UP,DLX,DLS, 
*IMAX , NMAX, NF 

1 IRST,NL AST, FACTO, TCOST,TPIV 
D  DIMENSIONS  FOR  COMMON 

DIMENSION  U(50> ,XK(60,60) ,S(60) ,X(60) ,Z(60) ,BK(60,60) , 
*B (60) , ZP (60) 

1, A (60, 60) , ALPHA (10) ,R( 60) ,XX(60) ,UP(60) 

UAV  =  0. 

DO  40  I  =  1 , 1  MAX 
COEF  =  0. 

DO  41  J  =  1 , NMA X 
KNUM  =  J/2 
JNUM  =  ( J*1 ) /2 
IF(JNUM.NE.KNUM)  GO  TO  43 
SIG  =  2. 

GO  TO  44 
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43  SIG  =  4. 

44  COEF  =  COEF  +  SIG*XK(I,J)*ZP(J)/(3®*DLS> 
4i  CONTINUE 

UP(I>  =  COEF 

uav  =  uav  +  (upm-um>»*2 

40  CONTINUE 

UAVG  =  SORT  (UAV/FLOAT ( XMA  X) ) 

90  RETURN 
END 


SUBROUTINE  REG2 
PURPOSE 

THIS  SUBROUTINE  INTEGRATES  XK(I,J)*XK(I, J1  OVER 
ZI-VALUES  TO  OBTAIN  NEW  MATRIX  BK(I,J) 

USAGE 

CALL  REG2 


SUBROUTINE  REG2 

COMMON/ZYT/U,XK,S,X,Z,BKsB,ZP,A,ALPHA,R»XX,UP,OLX,OLS, 

*IMAX,NMAX,NF 

1IRST, NL AST, FACTO, TCOST,TPIV 
DIMENSIONS  FOR  COMMON 

DIMENSION  U(BO) ,XK(60,60) , S (60) ,X(60),Z(60) ,BK (60,60) , 
*B(60) ,ZP(60> 

1, A (60, 60) , ALPHA (10) ,R(60> ,XX(60) ,UP(60) 

SIMPSON  RULE 

DO  5  I  =  1 , NM AX 
DO  5  J  =  1 , NM AX 
COEF 1  =  0. 

COEF  =  0. 

DO  20  K  =  1,IMAX 

IF(K.EQ.l.OR.K.EQ.IMAX)  GO  TO  21 
IF(K.EQ.2.0R.K.EQ.IMAX-1)  GO  TO  23 
KNUM  =  K/2 
JNUM  =  (K+l) t 2 
IF(JNUM.EQ.KNUM)  GO  TO  23 

22  SIG  =  2. 

GO  TO  24 

21  SIG  =  1. 

GO  TO  24 

23  SIG  =  4. 

24  A1  =  SIG*XK(K,I)*XK(K, J)/(3.*DLX) 

IF(I.GT.l)  GO  TO  7 
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A2  =  SIG*XK(K,J)*U(K) M3.*DLX) 

C0EF1  =  GOEFi  *  A2 

7  COEF  =  COEF  +  Ai 
20  CONTINUE 

IF(I.GT.l)  GO  TO  8 
B ( J)  =  COEFi 

8  BK(IfJ)  =  COEF/DLS 
5  CONTINUE 

RETURN 

END 

*««*****« ****»*»*»**» »**»*•*****»***»• **************** 


SUBROUTINE  REG3 
PURPOSE 

THIS  SUBROUTINE  INTRODUCES  THE  REGULARIZATION  TERMS 
IN  THE  MATRIX  BK(I,J>.  THE  FINAL  REGULARIZED 
MATRIX  IS  A(I,J) 

USAGE 

CALL  REG3 


SUBROUTINE  REG3 

COMMON/ZYT/U,XK,S,X,Z,BK,B,ZP,A,ALPHA,R,XX,UP,OLX,DLS, 

*IMAX,NMAX,NF 

i IRST  ,NL AST, FACTO,  TCOST  ,  TPIV 
C  DIMENSIONS  FOR  COMMON 

DIMENSION  U(60> ,XK (60,60) ,S(60> ,X(6Q) ,Z(60) ,BK (60,601  , 
*B(60> ,ZP(60) 

1, A (60,60) , ALPHA (10) ,R(60) ,XX(60) ,UP(60) 

DO  51  I  *  1 , NMAX 
DO  51  J  =  1, NMAX 
51  A (I , J)  =  0. 

DO  9  I  =  1,NMAX 
DO  9  J  =  1 , NMAX 
A  (I, J)  =  BK(I , J) 

9  CONTINUE 

DO  63  N  =  NFIRST,NLAST 
DO  60  I  =  1 , NMAX 
NUM  =  N  +  1 
DO  60  J  =  1 , NUM 
DO  60  K  =  1 , NUM 
NB  =  N/2 

LABEL1  =  I  -  N8  +  J  -  1 
LABEL2  *  I  •  NB  H  ■  1 

IF(LABELl»GT«NMAX»OR.L  ABEL  1. LT • 1)  GO  TO  60 
Al  =  CALC(N, J,DLS) 

A2  =  CALC(N,K,DLS) 
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50  IF(LABEL2.GT.NMAX.0R.LABEL2.LT.l)  GO  TO  59 
QQ  =  0.5*Al*A2*ALPHA(N)*(i.E  03*DLS)**N 
A(LABEL1,LA8EL2)  =  A(LABEL1,LABEL2»  ♦  QO 
GO  TO  60 

59  IF(L ABEL2. GT . NM AX)  GO  TO  58 
LABEL2  =  IABS  (LABEL2)  +  i 
GO  TO  50 

58  LABEL2  =  2*NMAX  -  LABEL2  +  1 
GO  TO  50 

60  CONTINUE 
63  CONTINUE 

RETURN 

END 

*»**»****»»**¥*****#» ******************* ****** **«»*¥** 


FUNCTION  CALC 
PURPOSE 

THIS  FUNCTION  SUBROUTINE  EVALUATES  THE 
COEFFICIENTS  (BINOMIAL),  ETC. 

CALLED  8T  REG3 

USAGE 

X=CALC(N,K,DLS) 


FUNCTION  CALC (N,K»DLS) 

LL  =  2*N 
L  =  K  -  1 
H  =  N  -  L 

IF(K.EQ.1.0R.K.E9.N+1)  GO  TO  10 

11  =  1 

12  =  1 
13  =  1 

DO  1  I  =  1 ,  L 

1  II  =  11*1 

DO  3  I  =  1 ,  N 
3  13  =  1*13 
DO  2  I  =  1 ,  H 

2  12  =  1*12 

XI  =  13/(11*12) 

X2  =  (-1. ) **K 
CALC  =  X1*X2 
GO  TO  99 

10  IF(K.EQ.l)  GO  TO  11 
XI  =  -1. 

CALC  =  XI 
GO  TO  99 
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il  XI  =  (-1,)**K 
CALC  =  Xi 
99  RETURN 
END 


SUBROUTINE  QUA01 
PURPOSE 

SUBROUTINES  QUAOi  AND  QUAD?  SOLVE  THE  QUADRATIC 
PROGRAMMING  PROBLEM,  FIND  THE  VALUE  OF  X  THAT 
MAXIMIZES 

A*X  -  1 12  X  *  B  X 

SUBJECT  TO 

C  X  .LE.  D 
X  .GE.  0 

USING  DANTZIG'S  MODIFIED  SIMPLEX  ALGORITHM 

(SEE  JOHN  C.  G.  BOOT,  QUADRATIC  PROGRAMMING,  RAND 

MCNALLV,  CHICAGO  1964,  PP.  186-196) 

QUAD1  DEFINES  THE  INITIAL  SIMPLEX  TABLEAU 

USAGE 
CALL  QUAD1 

( A, B,C,D,N,K,ROWS, COLS, ABCD, RESULT, ZER01, BASIS) 
DESCRIPTION  OF  PARAMETERS 

A  -  INPUT  VECTOR  OF  LENGTH  N  THAT  DEFINES 

THE  LINEAR  PART  OF  THE  OBJECTIVE  FUNCTION 
B  -  INPUT  MATRIX  ( N,N)  THAT  DEFINES  THE 

QUADRATIC  PART  OF  THE  OBJECTIVE  FUNCTION 
C  -  INPUT  MATRIX  (K,N)  THAT  DEFINES  THE  LEFT 

HAND  PART  OF  THE  CONSTRAINTS 
D  -  VECTOR  OF  LENGTH  K  THAT  DEFINES  THE  RIGHT 

HAND  SIDE  OF  THE  CONSTRAINTS 
N  -  NUMBER  OF  ELEMENTS  OF  X 

K  -  NUMBER  OF  CONSTRAINTS 

ROWS  -  N+K,  THE  NUMBER  OF  ROWS  IN  THE  INITIAL 

SIMPLEX  TABLEAU 

COLS  -  2*R0WS+i,  THE  NUMBER  OF  COLUMNS  IN  THE 

INITIAL  SIMPLEX  TABLEAU 

ABC A  -  THE  INITIAL  SIMPLEX  TABLEAU,  MATRIX  OF 

SIZE  (ROWS, COLS) 

RESULT  -  VECTOR  OF  LENGTH  ROWS,  THAT  CONTAINS  THE 
RESULTS  OF  THE  QUADRATIC  PROGRAMMING 
PROBLEM 

BASIS  -  VECTOR  OF  LENGTH  ROWS,  CONTAINING  THE 
LOCATIONS  OF  THE  BASIS  VECTORS 
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THE  BASIS  VECTORS 


SUBROUTINE  OUAOKRESULT, BASIS) 

INTEGER  ROHS, COLS, ZEROl, BASIS 

COMMON/ZYT/U,XK,S,X,Z,BK,B,ZP,A,ALPHA,R,XX,UP»DLX,OLS, 

*IMAX,NMAX,NF 

1IRST,NLAST,FAGT0,TC0ST,TPIV 
C  DIMENSION  FOR  COMMON 

DIMENSION  U(60) ,XK(60» 60) ,S(60) ,X(60) ,Z(60) ,BK(60,60) , 
*B(60) ,ZP(60) 

1, A (60,60) , ALPHA (10) , R(60) ,XX(60) ,UP(60) 

COMMON/ZXT/ABCO 

DIMENSION  C(2,60> ,D(2) , ABCD(6 0, 121) , RESULT (60) , ZEROl (1 
*21) , BASIS(60 
1) ,IROH (121) 

LOGICAL  NOPIVT 
COMMON/QUAD2C/NOPIVT 
N  =  NMAX 
K  =  1 

ROHS  =  NMAX  ♦  K 
COLS  =  2* (NMAX  ♦  Kl  ♦  1 
DO  57  I  *  1,NMAX 
C(1,I)  =  FACTO/OLS 
57  CONTINUE 

0(1)  =  FACTO 

C  N  VARIBLES  K  CONSTRAINTS 
DO  1  1=1, ROHS 
DO  1  J=l,GOLS 

1  ABCD(I,J)=0.0 
DO  2  1=1, N 

DO  2  J=i»N 

2  ABCD (I , J)  =  -A (I, J) 

00  3  Ki=i, K 
I=N+K1 

ABCD (1,1) =i»  0 
DO  3  J=l, N 

3  ABCD (I , J) =G (Kl,  J) 

DO  4  1=1, N 
J=N+K+I 

4  ABCD (I , J) =1.0 
DO  5  K1=1,K 
J=2*N+K+K1 

DO  5  1=1, N 

5  ABCD (I ,J)=-C(K1,I) 

J=COLS 

DO  6  1=1, N 

5  ABCD (I , J)  =  -B(I) 

DO  7  K1=1,K 
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I=N+K1 

I  ABCD(I,J)=D(Ki> 

DO  11  1=1, ROMS 
DO  12  J=1 , COLS 

12  IROW(J)=ABCO(I, J) 

II  CONTINUE 

CALL  QUAD2 (ROMS, COLS, N,K, RESULT, ZEROl, BASIS, TPIV,TC0ST 
*) 

RETURN 
END 


SUBROUTINE  QUAD2 

PURPOSE 

(SEE  QUAD1) 

USAGE 

CALL  0UA02 ( ABCD, ROWS, COLS, N,K, RESULT, ZEROl, BASIS) 

DESCRIPTION  OF  PARAMETERS 
(SEE  QUAOl) 


SUBROUTINE  0UAD2 (ROMS, COLS ,N,K, RESULT, ZEROl .BASIS, TPIV 
* , TCOST) 

INTEGER  COLS, ZEROl, ROMS, BASIS, COLN,PIVROM, PC, PR 
COMMON/ZXT/ABCD 

DIMENSION  ABCD (60 , 121) , RESULT (60) , ZEROl (121) ,BASIS(60) 
REAL  NUM, MULT 

LOGICAL  NOPIVT 
C OMM ON /QUAD2C /NOPIVT 
C  CLEAR  ZEROl  VECTOR 
DO  2  1=1, COLS 

2  ZEROl ( I) =0 

C  INSERT  (N+K)  ONES  INTO  ZEROl (N+1) 

DO  3  1=1, ROWS 
J=N+I 

3  ZEROl ( J) =1 

C  LOAD  N  COLUMN  NUMBERS  FROM  VARIABLE  V(l)  INTO  BASIS(i) 

DO  4  1=1, N 

4  BASIS(I)=ROWS+I 

C  LOAD  K  COLUMN  NUMBERS  FROM  VARIABLE  L(N+1)  INTO 
C  BASIS (N+l) 

DO  9  1=1, K 
J=N+I 

9  B ASIS ( J) = J 

C  ASSUME  A  NON  STANDARD  TABLE 
1?  NONSTD=0 
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3  LOOK  AT  ZEROi  VECTOR  AND  DETERMINE  FOR  A  NONSTANDARD 
S  TABLE . 

C  1.  PC=COLUMN  NUMBER  OF  THE  MISSING  V  VARIABLE 

C  2.  IV=V  COLUMN  NUMBER  OF  THE  BASIC  PAIR 

C  3.  IF  CONDITION  1  ANO  2  ARE  PRESENT  SET  NONSTO^l 

00  5  1=1, ROWS 
J=I*ROWS 

Il=ZEROi (I) + ZEROI (  J) 

IF  <I1-1)  6,5,7 
S  PC=J 

N0NST0=1 
GO  TO  5 

7  IV=J 
N0NSTD=1 

5  CONTINUE 

C  IS  THIS  A  NON  STANDARD  TABLE 
IF  (NONSTD.EQ.i)  GO  TO  8 

C  SCAN  THE  BASIS  FOR  IV, COLUMN  NUMBER  OF  THE  LARGEST 
0  NEGATIVE  V(I)  AND  DETERMINE  PC=COLUMN  NUMBER  OF  L(I)  TO 
C  BE  ADDED  TO  THE  BASIS 
VNEG=0 • 0 
DO  10  1=1, ROWS 
COLN=BASIS (I) 

IF  (COLN.LE.ROWS)  GO  TO  10 
T1=ABC0(I,G0LS) 

IF  (Tl.GE.VNEG)  GO  TO  10 
VNEG=T 1 
IV=COLN 
PC=COLN-ROWS 
10  CONTINUE 

C  LOOK  AT  THE  V(IV)  RATIO  AND  ALL  JCI)  RATIOS  AND 
C  DETERMINE  THE  VARIABLE  HAVING  THE  SMALLEST  NON  NEGATIVE 
C  VALUE.  THIS  COLUMN  NUMBER  IS  PR  (THE  PIVOT  ROW),  THE 
C  VARIABLE  TO  BE  REMOVED 

8  RATI0=1. 0E37 

NOPIVT= .TRUE. 

DO  11  1=1, ROWS 
COLN=BASIS(I) 

DEN= ABCD( I , PC) 

NUM=  ABCD( I , COLS) 

IF ( A3SCDEN)  .LT.  TPIV)GO  TO  11 
IF  (COLN.LE.ROWS)  GO  TO  13 
IF  (COLN.NE. IV)  GO  TO  11 
13  T 1=NUM/DEN 

NOPIVT=. FALSE. 

IF  (Tl.LE.O.O)  GO  TO  11 
IF  (Tl.GE. RATIO)  GO  TO  11 
PR=COLN 
PIVROW=I 
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RATIO=Ti 

li  CONTINUE 

0  ADD  ANO  DELETE  THE  PROPER  VARIABLES  FROM  THE  BASIS  AND 

0  ZEROi  VECTORS 
ZEROi (PC) =1 
ZEROi (PR) =0 
BASIS (PIVROW) =PC 
PR=PIVROW 

0  NORMALIZE  THE  PIVOT  ROW  BY  THE  PIVOT  ELEMENT 
DEN= ABCO (PR, PC) 

DO  14  J=i,COLS 

14  ABCO (PR*J>=A8CD(PR,J>  /  DEN 

C  ZERO  OUT  THE  REMAINING  ELEMENTS  OF  THE  PIVOT  COLUMN 
DO  18  1=1, ROWS 
IF  (I.EQ.PR)  GO  TO  18 
MULT=-ABCO (I, PC) 

IF  (MULT.EO.O.O)  GO  TO  18 
DO  15  J=i , COLS 

15  ABCO (I»J)=ABCD(I,J)+MULT*ABCD(PR,J) 

18  CONTINUE 

c  are  any  of  the  basic  values  still  negative 

DO  16  1=1, ROWS 

IF(ABCD (I, COLS)  .LT.  -ABS (TCOST) ) GO  TO  17 

16  CONTINUE 

C  TRANSFER  THE  LAST  COLUMN  TO  THE  SOLUTION  VECTOR 

DO  1  1=1, ROWS 

1  RESULT (I)=ABCD( I, COLS) 

RETURN 

END 
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THIS  IS  FOR  A  SYMMETRICAL  UNIMODAL  MOL.  HEIGHT  DIST.  USIN 
3  REQUAD 
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The  equation  relating  molecular  weight  distribution  of  a  polymer  to  the  experi¬ 
mental  function  of  concentration  appearing  in  quilibrium  sedimentation  with  the 
ultracentrifuge  is  nonsol vable  because  it  is  an  Improperly  Posed  Problem  in  the 
Hadamard  sense.  For  a  simple  distribution  this  equation  has  been  solved  by  applying 
a  method  of  regularization.  To  solve  a  nonsymmetrical  bimodal  and  a  trimodal  dis¬ 
tribution,  the  technique  of  regularization  had  to  be  incorporated  into  a  linear  pro¬ 
gramming.  In  the  current  work  the  regularization  technique  has  been  incorporated  into 
quadratic  programming.  This  new  combined  method  proved  to  be  more  adequate  to  solve, 
also  more  complex  distributions  such  as  tri-,  tetra-,  and  pentamodal .  In  addition, 
this  technique  is  cheaper,  because  it  requires  less  computer  time  than  the  regular¬ 
ization  incorporated  into  linear  programming. 
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